Load packages

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.1.0     ✓ dplyr   1.0.5
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

Check working directory

getwd() 
## [1] "/Users/jillashey/Desktop/PutnamLab/Repositories/Astrangia_repo/Astrangia_repo"

Load data

# Load data
daily <- read.csv("data/DailyMeasurements.csv", header = T)
daily <- daily[-13,] # remove tank 14 because it cracked 
daily$Timepoint <- as.Date(daily$Timepoint, "%m/%d/%Y")
daily$Tank <- as.factor(daily$Tank)
daily <- subset(daily, Treatment=="Heat" | Treatment=="Ambient")

Discrete pH calculations from Tris calibrations.

path <-("data/pH_Tris") #set path to calibration files
file.names<-list.files(path = path, pattern = "csv$") #list all the file names in the folder to get only get the csv files
pH.cals <- data.frame(matrix(NA, nrow=length(file.names), ncol=3, dimnames=list(file.names,c("Date", "Intercept", "Slope")))) #generate a 3 column dataframe with specific column names

for(i in 1:length(file.names)) { # for every file in list start at the first and run this following function
  Calib.Data <-read.table(file.path(path,file.names[i]), header=TRUE, sep=",", na.string="NA", as.is=TRUE) #reads in the data files
  file.names[i]
  model <-lm(mVTris ~ TTris, data=Calib.Data) #runs a linear regression of mV as a function of temperature
  coe <- coef(model) #extracts the coeffecients
  summary(model)$r.squared #extracts the r squared
  plot(Calib.Data$mVTris, Calib.Data$TTris) #plots the regression data
  pH.cals[i,2:3] <- coe #inserts coefficients in the dataframe
  pH.cals[i,1] <- substr(file.names[i],1,8) #stores the file name in the Date column
}

colnames(pH.cals) <- c("Calib.Date",  "Intercept",  "Slope") #rename columns
pH.cals #view data
##              Calib.Date Intercept      Slope
## 20210122.csv   20210122 -86.73848  1.0218902
## 20210126.csv   20210126 -86.23930  0.9322418
## 20210203.csv   20210203 -88.14064  1.0017943
## 20210205.csv   20210205 -88.23311  0.8416322
## 20210217.csv   20210217 -87.53924  1.1237563
## 20210221.csv   20210221 -85.87848  0.8968835
## 20210223.csv   20210223 -86.36221  0.9370372
## 20210225.csv   20210225 -85.86760  0.9234509
## 20210227.csv   20210227 -87.60427  1.0150785
## 20210308.csv   20210308 -89.42691  0.9977500
## 20210317.csv   20210317 -86.21007  0.9799588
## 20210325.csv   20210325 -87.82550  1.0238389
## 20210409.csv   20210409 -88.31599  1.1153461
## 20210411.csv   20210411 -53.62103 -0.1385311
## 20210420.csv   20210420 -68.46091  0.2421390
## 20210524.csv   20210524 -96.95526  1.1877897
## 20210527.csv   20210527 -98.21886  1.2378231
## 20210531.csv   20210531 -88.25076  0.8102227
#constants for use in pH calculation 
R <- 8.31447215 #gas constant in J mol-1 K-1 
F <- 96485.339924 #Faraday constant in coulombs mol-1

#merge with Seawater chemistry file
SW.chem <- merge(daily, pH.cals, by="Calib.Date")

Calculate total pH.

mvTris <- SW.chem$Temperature*SW.chem$Slope+SW.chem$Intercept #calculate the mV of the tris standard using the temperature mv relationships in the measured standard curves 
STris<-35 #salinity of the Tris
phTris<- (11911.08-18.2499*STris-0.039336*STris^2)*(1/(SW.chem$Temperature+273.15))-366.27059+ 0.53993607*STris+0.00016329*STris^2+(64.52243-0.084041*STris)*log(SW.chem$Temperature+273.15)-0.11149858*(SW.chem$Temperature+273.15) #calculate the pH of the tris (Dickson A. G., Sabine C. L. and Christian J. R., SOP 6a)
SW.chem$pH.Total<-phTris+(mvTris/1000-SW.chem$pH.MV/1000)/(R*(SW.chem$Temperature+273.15)*log(10)/F) #calculate the pH on the total scale (Dickson A. G., Sabine C. L. and Christian J. R., SOP 6a)

Plot daily measurements

# By Treatment
temp.trt = ggplot(SW.chem, aes(x=Treatment, y=Temperature.C)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("Temperature°C") +
  theme(axis.text.x = element_text(angle = 90))
temp.trt
## Warning: Removed 9 rows containing non-finite values (stat_boxplot).

sal.trt = ggplot(SW.chem, aes(x=Treatment, y=Salinity.psu)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("Salinity.psu") +
    theme(axis.text.x = element_text(angle = 90))
sal.trt
## Warning: Removed 21 rows containing non-finite values (stat_boxplot).

pH.trt = ggplot(SW.chem, aes(x=Treatment, y=pH.Total)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("pH Total Scale") +
  theme(axis.text.x = element_text(angle = 90))
pH.trt
## Warning: Removed 22 rows containing non-finite values (stat_boxplot).

light.trt = ggplot(SW.chem, aes(x=Treatment, y=Light)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("Light µmol photons m^2 s^-1") +
    theme(axis.text.x = element_text(angle = 90))
light.trt
## Warning: Removed 625 rows containing non-finite values (stat_boxplot).

flow.trt = ggplot(SW.chem, aes(x=Treatment, y=Flow.mL.sec)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("Flow mL/sec") +
    theme(axis.text.x = element_text(angle = 90))
flow.trt
## Warning: Removed 1040 rows containing non-finite values (stat_boxplot).

plot.trt <- grid.arrange(temp.trt, sal.trt, pH.trt, light.trt, flow.trt, ncol=3, nrow = 2, clip="off")
## Warning: Removed 9 rows containing non-finite values (stat_boxplot).
## Warning: Removed 21 rows containing non-finite values (stat_boxplot).
## Warning: Removed 22 rows containing non-finite values (stat_boxplot).
## Warning: Removed 625 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1040 rows containing non-finite values (stat_boxplot).

ggsave("output/Daily_Measurements.By.Treatment.pdf", plot.trt, width = 21, height = 21, units = c("in"))

# By Tank
temp.tank = ggplot(SW.chem, aes(x=Tank, y=Temperature.C)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("Temperature°C") +
    theme(axis.text.x = element_text(angle = 90))
temp.tank
## Warning: Removed 9 rows containing non-finite values (stat_boxplot).

sal.tank = ggplot(SW.chem, aes(x=Tank, y=Salinity.psu)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("Salinity.psu") +
    theme(axis.text.x = element_text(angle = 90))
sal.tank
## Warning: Removed 21 rows containing non-finite values (stat_boxplot).

pH.tank = ggplot(SW.chem, aes(x=Tank, y=pH.Total)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("pH Total Scale") +
    theme(axis.text.x = element_text(angle = 90))
pH.tank
## Warning: Removed 22 rows containing non-finite values (stat_boxplot).

light.tank = ggplot(SW.chem, aes(x=Tank, y=Light)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("Light µmol photons m^2 s^-1") +
    theme(axis.text.x = element_text(angle = 90))
light.tank
## Warning: Removed 625 rows containing non-finite values (stat_boxplot).

flow.tank = ggplot(SW.chem, aes(x=Tank, y=Flow.mL.sec)) +
  geom_boxplot(aes(color = Treatment), size = 1) +
  ylab("Flow mL/sec") +
    theme(axis.text.x = element_text(angle = 90))
flow.tank
## Warning: Removed 1040 rows containing non-finite values (stat_boxplot).

plot.tank <- grid.arrange(temp.tank, sal.tank, pH.tank, light.tank, flow.tank, ncol=3, nrow = 2, clip="off")
## Warning: Removed 9 rows containing non-finite values (stat_boxplot).
## Warning: Removed 21 rows containing non-finite values (stat_boxplot).
## Warning: Removed 22 rows containing non-finite values (stat_boxplot).
## Warning: Removed 625 rows containing non-finite values (stat_boxplot).
## Warning: Removed 1040 rows containing non-finite values (stat_boxplot).

ggsave("output/Daily_Measurements.By.Tank.pdf", plot.tank, width = 21, height = 21, units = c("in"))

test for differences between tanks and treatments

# By Treatment
mod.temp <- aov(Temperature.C ~ Treatment, data = SW.chem)
mod.temp
## Call:
##    aov(formula = Temperature.C ~ Treatment, data = SW.chem)
## 
## Terms:
##                 Treatment Residuals
## Sum of Squares   3375.159  8500.317
## Deg. of Freedom         1      1109
## 
## Residual standard error: 2.768546
## Estimated effects may be unbalanced
## 9 observations deleted due to missingness
summary(mod.temp)
##               Df Sum Sq Mean Sq F value Pr(>F)    
## Treatment      1   3375    3375   440.3 <2e-16 ***
## Residuals   1109   8500       8                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 9 observations deleted due to missingness
hist(mod.temp$residuals)

mod.sal <- aov(Salinity.psu ~ Treatment, data = SW.chem)
mod.sal
## Call:
##    aov(formula = Salinity.psu ~ Treatment, data = SW.chem)
## 
## Terms:
##                 Treatment Residuals
## Sum of Squares    18.7864 1163.1030
## Deg. of Freedom         1      1097
## 
## Residual standard error: 1.029688
## Estimated effects may be unbalanced
## 21 observations deleted due to missingness
summary(mod.sal)
##               Df Sum Sq Mean Sq F value   Pr(>F)    
## Treatment      1   18.8   18.79   17.72 2.77e-05 ***
## Residuals   1097 1163.1    1.06                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 21 observations deleted due to missingness
hist(mod.sal$residuals)

mod.pH <- aov(pH.Total ~ Treatment, data = SW.chem)
mod.pH
## Call:
##    aov(formula = pH.Total ~ Treatment, data = SW.chem)
## 
## Terms:
##                 Treatment Residuals
## Sum of Squares    2.57166  32.70446
## Deg. of Freedom         1      1096
## 
## Residual standard error: 0.1727421
## Estimated effects may be unbalanced
## 22 observations deleted due to missingness
summary(mod.pH)
##               Df Sum Sq Mean Sq F value Pr(>F)    
## Treatment      1   2.57  2.5717   86.18 <2e-16 ***
## Residuals   1096  32.70  0.0298                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 22 observations deleted due to missingness
hist(mod.pH$residuals)

mod.light <- aov(Light ~ Treatment, data = SW.chem)
mod.light
## Call:
##    aov(formula = Light ~ Treatment, data = SW.chem)
## 
## Terms:
##                 Treatment Residuals
## Sum of Squares      90.49 124593.05
## Deg. of Freedom         1       493
## 
## Residual standard error: 15.8973
## Estimated effects may be unbalanced
## 625 observations deleted due to missingness
summary(mod.light)
##              Df Sum Sq Mean Sq F value Pr(>F)
## Treatment     1     90   90.49   0.358   0.55
## Residuals   493 124593  252.72               
## 625 observations deleted due to missingness
hist(mod.light$residuals)

mod.flow <- aov(Flow.mL.sec ~ Treatment, data = SW.chem)
mod.flow
## Call:
##    aov(formula = Flow.mL.sec ~ Treatment, data = SW.chem)
## 
## Terms:
##                 Treatment Residuals
## Sum of Squares     3167.9  679148.5
## Deg. of Freedom         1        78
## 
## Residual standard error: 93.31148
## Estimated effects may be unbalanced
## 1040 observations deleted due to missingness
summary(mod.flow)
##             Df Sum Sq Mean Sq F value Pr(>F)
## Treatment    1   3168    3168   0.364  0.548
## Residuals   78 679148    8707               
## 1040 observations deleted due to missingness
hist(mod.flow$residuals)

# By Tank
mod.temp <- aov(Temperature.C ~ Tank, data = SW.chem)
mod.temp
## Call:
##    aov(formula = Temperature.C ~ Tank, data = SW.chem)
## 
## Terms:
##                     Tank Residuals
## Sum of Squares  3384.931  8490.545
## Deg. of Freedom       15      1095
## 
## Residual standard error: 2.784587
## Estimated effects may be unbalanced
## 9 observations deleted due to missingness
summary(mod.temp)
##               Df Sum Sq Mean Sq F value Pr(>F)    
## Tank          15   3385  225.66    29.1 <2e-16 ***
## Residuals   1095   8491    7.75                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 9 observations deleted due to missingness
hist(mod.temp$residuals)

mod.sal <- aov(Salinity.psu ~ Tank, data = SW.chem)
mod.sal
## Call:
##    aov(formula = Salinity.psu ~ Tank, data = SW.chem)
## 
## Terms:
##                      Tank Residuals
## Sum of Squares    23.8552 1158.0342
## Deg. of Freedom        15      1083
## 
## Residual standard error: 1.034062
## Estimated effects may be unbalanced
## 21 observations deleted due to missingness
summary(mod.sal)
##               Df Sum Sq Mean Sq F value Pr(>F)
## Tank          15   23.9   1.590   1.487  0.102
## Residuals   1083 1158.0   1.069               
## 21 observations deleted due to missingness
hist(mod.sal$residuals)

mod.pH <- aov(pH.Total ~ Tank, data = SW.chem)
mod.pH
## Call:
##    aov(formula = pH.Total ~ Tank, data = SW.chem)
## 
## Terms:
##                     Tank Residuals
## Sum of Squares   2.71815  32.55798
## Deg. of Freedom       15      1082
## 
## Residual standard error: 0.1734663
## Estimated effects may be unbalanced
## 22 observations deleted due to missingness
summary(mod.pH)
##               Df Sum Sq Mean Sq F value   Pr(>F)    
## Tank          15   2.72 0.18121   6.022 3.19e-12 ***
## Residuals   1082  32.56 0.03009                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 22 observations deleted due to missingness
hist(mod.pH$residuals)

mod.light <- aov(Light ~ Tank, data = SW.chem)
mod.light
## Call:
##    aov(formula = Light ~ Tank, data = SW.chem)
## 
## Terms:
##                      Tank Residuals
## Sum of Squares    7884.24 116799.31
## Deg. of Freedom        15       479
## 
## Residual standard error: 15.61537
## Estimated effects may be unbalanced
## 625 observations deleted due to missingness
summary(mod.light)
##              Df Sum Sq Mean Sq F value  Pr(>F)   
## Tank         15   7884   525.6   2.156 0.00705 **
## Residuals   479 116799   243.8                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 625 observations deleted due to missingness
hist(mod.light$residuals)

mod.flow <- aov(Flow.mL.sec ~ Tank, data = SW.chem)
mod.flow
## Call:
##    aov(formula = Flow.mL.sec ~ Tank, data = SW.chem)
## 
## Terms:
##                     Tank Residuals
## Sum of Squares   51655.4  630661.0
## Deg. of Freedom       15        64
## 
## Residual standard error: 99.26771
## Estimated effects may be unbalanced
## 1040 observations deleted due to missingness
summary(mod.flow)
##             Df Sum Sq Mean Sq F value Pr(>F)
## Tank        15  51655    3444   0.349  0.987
## Residuals   64 630661    9854               
## 1040 observations deleted due to missingness
hist(mod.flow$residuals)

Plot by date

# Plot and save
# By Treatment
timepoint.trt <- ggplot(SW.chem, aes(x=Timepoint, y=Temperature.C, group=Treatment)) +
  geom_line(aes(color = Treatment), size = 1) 

# By Tank
timepoint.tank <- ggplot(SW.chem, aes(x=Timepoint, y=Temperature.C)) +
  geom_line(aes(color = Tank), size = 1) +
  scale_fill_brewer("Dark2") 

plot.timepoint <- grid.arrange(timepoint.tank, timepoint.trt, ncol=1, nrow = 2, clip="off")

ggsave("output/Daily_Measurements.By.Date.pdf", plot.tank, width = 21, height = 25, units = c("in"))